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Abstract 

We explore the effect of higher order operators in the non-relativistic formulation of 
QCD (NRQCD). We calculated masses in the 56-spectrum using quenched gauge con¬ 
figurations at /3 = 6.0 and two different NRQCD actions which have been corrected to 
order mv^ and mv^. The two-point functions are calculated in a gauge invariant fashion. 

We find the general structure of the spectrum to be the same in the two cases. Using 
the IP — l^S'i splitting we determine the inverse lattice spacings to be 2.44(4) GeV and 
2.44(5) GeV for the mv'^ and mv^ actions, respectively. We do observe shifts in the spin 
splittings. The hyperfine splitting is reduced by approximately 4 MeV, while the fine 
splitting is down by up to 10 MeV, albeit with large statistical errors. 


1 Introduction 

Heavy quark systems have long been studied and much experimental data has been accumu¬ 
lated which will eventually lead to tight constraints on the parameters of the Standard Model. 
To this end it is important to have accurate non-perturbative predictions from QCD which 
can be compared with experiment. One such method is Lattice QCD and it is necessary to 
understand how heavy quark systems can be treated in this framework. 

A nonrelativistic approximation to QCD, NRQCD, was proposed to go beyond the 
static approximation in lattice calculations. NRQCD has been remarkably successful in re¬ 
producing the spectrum of heavy quark systems |Q, ^ owing to the fact that the quarks 

within such hadrons move with velocities v such that . Furthermore, NRQCD has the 

virtue of retaining, at least approximatively, the quark dynamics and can be considered an 
effective theory. Computationally, such an approach is much faster to solve than the inversion 
problem in fully relativistic QCD. A systematic improvement program has been developed to 
match lattice NRQCD to continuum QCD 1,0- As an effective field theory the predictive 
power of NRQCD relies on the control of higher dimensional operators. 

In this paper we calculate the bb spectrum comparing two different NRQCD-actions with 
accuracies 0{mv^) and 0{mv^). In section 2 we give the details of our evolution equation 
and the gauge invariant operators used. In section 3 we state the results and convert them 
into dimensionful units. In the concluding section we compare our results with other work in 
this area and experimental data. 
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2 NRQCD and Operators 


In NRQCD the inversion problem of the fermion matrix is an initial value problem. We use 
the evolution equation along the Euclidean time-direction defined by: 
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G{x,ty + l;y) 
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Here S(x,y) is the source at the first timeslice (t = ty). We have S(x,y) = S^^^(x,y) for 
a single quark source at the origin, y, but we also propagate extended objects with the 
same evolution equation. The operator Hq is the leading kinetic term and 5H contains the 
relativistic corrections and spin corrections. The last two terms in 6H are present to correct 
for lattice spacing errors in temporal and spatial derivatives. For the derivatives we give the 
following definitions which are consistent with those of |P, ^]; 
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For our calculation with accuracy 0{mv^) we set 04,05 and cq to zero and replace A by A. 
In this case we calculate the helds E and B from the clover field Fyi, in the standard fashion 
1^]. To determine the spin splittings with accuracy 0{mv^) we retain the above coefficients 
and the improved operator A. Correspondingly, we replace standard discretized gauge field 
term by 

Fyu = ^Fy^-^(Uf,{x)Fy^{x + fj,)Ul{x) + U-y{x)Fy^{x-y)Ul^{x)-{fi^iy)^ , (4) 
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when this accuracy is needed for E and B in Equation |^. 

All the gauge links are tadpole improved with uq = (0|^Trf7^i,|0)^/^ as suggested in [^. It 
has been demonstrated that tadpole improvement is crucial to reproduce the spin splittings 
accurately, which are otherwise underestimated. We assume that radiative corrections are 
sufficiently accounted for by tadpole improvement and we therefore set all the renormalisation 
coefficients, Cj, to 1. 

To extract masses we calculate two-point functions of operators with the appropriate 
quantum numbers. In a non-relativistic setting gauge-invariant meson operators can be con¬ 
structed from the two-spinors and 'l'(y) (which represent the antiquark field and quark 

field) and a Wilson line, W (x, y). Such operators have the following behaviour under rotation, 
parity and charge conjugation 

x\x)W{x,y)^{y) ^ x\Rx)W^{R x, Ry)^{Ry) , 

x\x)W{x,y)^{y) -x^(xp)W'^(xp, r/p)T(yp) , 

xHx)W{x,y)'if{y) -^'^{x)W^{x,y)x*{y) ■ (5) 


The construction of spatially extended operators on the lattice has been reported elsewhere 
1^]. Here we adopt the following notation: 

L”(x) = Ui{x)Ui{x + i)... Ui{x + {n - l)i) = L'^^{x + in) , 

A”T(x) = Lf (x)T(x-b m) — L”j(x)'I'(x — m) . (6) 


With the generalised link variables of equation and the transformation properties of equation 
^ we can construct operators with definite , where J labels the irreducible representations 
of the octahedral group (J = Ai, A 2 , Ti, T 2 , E). 

In order to increase the overlap between the meson and the ground states we create on 
the lattice we use extended versions of the derivative as defined in equation The operators 
we used are listed in Table |||. The meson correlator is written as a Monte Carlo average over 
all configurations 


C--(x,y) = (tr Gt(x,y)G--(x,y) ) 


(7) 


where tr denotes contraction over all internal degrees of freedom and G”™' is the smeared 
propagator defined by 


G--(x,y) ^ ^ 0-(x,zi)G(zi,Z2)0™t(^2,y) . (8) 

21,22 

Here (n, m) stands for the radii at (sink, source) and O” is an operator as defined in Table ^ 
but with extended symmetric derivative A"’. For the smeared propagator we solve equation 
1^ with S'(x,y) = 0”^^(x, y) and multiply with O” at the sink. We fix the origin at some 
(arbitrary) lattice point, y, and sum over all spatial x so as to project out the zero momentum 
mode. In addition we sum over all polarisations to increase the statistics. 

3 Simulation and Results 

We use quenched gauge field configurations at /3 = 6.0 on a 16^ x 48 lattice generated with the 
standard Wilson action at the EPCC in Edinburgh. The propagators were calculated at the 
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HPCF in Cambridge and at the Hitachi Europe Ltd. High Performance Computing Center in 
Maidenhead. As the operators chosen are gauge invariant, gauge fixing is not necessary. We 
used the extended operators defined in section 2 with radii 1,2 and 3 in both the source and 
sink thus giving rise to a 3 x 3 matrix of correlators for each set of quantum numbers. The 
calculation were done with a bare quark mass of ami, = 1.71 which is the same as that used 
in 1^. We used n = 2 as the stability parameter in Equation The tadpole improvement 
coefficient appropriate to j3 = 6.0 is uq = 0.878 as in Q. We calculated quark propagators for 
both actions from 8 points (i, j, k) where i, j, A: is 0 or 7 on the first timeslice for each of the 499 
configurations. It has been noted previously ||5| that such an arrangement gives independent 
and uncorrelated measurements due to the small size of the 65-system on our lattice. Thus 
we have a total of 35928 correlators for each meson we considered (^iS'o,^5i, ^Pi, ^Pb,i, 2 )- We 
also average over all polarisations and all spin components. 

We fit the correlators to the multi-exponential form 

iiflt 

= . ( 9 ) 

i=l 

Here a denotes a meson state with certain quantum numbers, (n, m) the different radii at the 
(sink, source) and t the Euclidean time. We chose ngt to be 1,2 or 3 and restrict ourselves 
to some conveniently chosen subsets of all the data, e.g. a window [tminjtmax]- The basic 
method for obtaining the parameters a”™ and M“ (which we generically refer to as {6;}) 
is a correlated ht using the whole covariance matrix. The statistical nature of the data 
causes fluctuations from the central values of the parameters, which are determined from a 
minimal x^{b) ht. We estimate the covariance matrix for the parameters from the inverse of 
/{dhkdbi). The goodness of the fit is quantified by the Q-value, defined in 0 - We 
require an acceptable fit to have Q > 0.1 and that the results remain consistent as we alter 
our fitting prescription slightly (e.g. varying tmin or the number of eigenvalues retained when 
inverting the covariance matrix of the data). To illustrate the method we display a fitted 
mass plotted against tmin in Eigure [T|. Our dimensionless results for the non-relativistic rest 
energies are shown in Table |2|. 

Eitting excitation energies relative to the ground state makes it difficult to extract the 
splitting and the associated error of the P-state. We rewrite Equation ^ as 

(^"'"(t) = , (10) 

i=\ 

where we can now fit the amplitudes a^, the ground state mass and the splittings, 5Mf-, 
with respect to Ml. Our results are tabulated in Table ^ 

Erom the IP — l^Pi splitting we estimate the lattice spacing at /3 = 6.0 to be = 

2.44(4) GeV for the improved action and = 2.44(5) GeV for the unimproved action. 

The kinetic masses were determined using 80 configurations for both actions. We show an 
example of the dispersion relation in Figure || and the results are shown in Table ^ 
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4 Conclusions 


Despite the success of NRQCD, a source of concern has been the effect of higher dimensional 
operators in C^jiqcd- To test this issue, we calculated the T spectrum in a gauge invariant 
fashion to accuracy 0{mv^). We then retained the spin corrections in Equation || up to 
0{mv^) and compared the results with each other. We also checked higher order contributions 
to the derivative in the term proportional to C 3 and found their effect negligible. Our main 
results are summarised in Figures and ^ In Table ^ we compare the results for the two 
different evolution equations with each other and convert the dimensionless numbers into 
physical units. We set the scale using the IP — splitting and find = 2.44(4) GeV 
and 2.44(5) GeV, respectively. Our results for a lower order evolution equation match those 
from an earlier calculation |P, |l^, U] where the configurations were fixed in the Goulomb 


gauge. 

These comparisons demonstrate the consistency of NRQCD. The gross structure of the 
T-spectrum calculated is unchanged when we added the spin correction terms. For example, 
the ratio Rsp = (2^S'i — 1^5i)/(l.P — 1^5i) = 1.3(1) agrees for both actions and with the 
experimental value, 1.28, within statistical errors. 

The bare mass parameter, airib, was chosen to be 1.71 to match a choice in [^. With our 
lattice spacing this corresponds to mb = 4.17(7) GeV and gives a kinetic mass of 9.4(3) GeV 
(for the mv^ action) and 9.7(3) GeV (for the mv'^ action) which should be compared to the 
experimental value of 9.46 GeV. The differences in the kinetic mass due to the extra terms 
is of the order of 1-2 standard deviations. We note also that the x^/dof resulting from the 
kinetic mass fits for the mv^ action is much smaller than that for the action. 

However, differences can be seen for the spin splittings. The introduction of the extra 
terms reduces the hyperfine splitting, — ^5o, by approximately 4 MeV. This size of a shift 
is what one would expect from power counting arguments. 

We also observe a reduction in the fine structure. For example, by adding the spin 
corrections the ^i -2 ~ ^Pi splitting is reduced from 19(3) MeV to 11(3) MeV, compared with 
the experimental value of 21 GeV. Such a behaviour has also been reported for the spectrum 
calculations of charmonium on much coarser lattices |]l^. Here we find a similar tendency 
and further analysis will be needed to decide whether this feature persists in unquenched 
calculations. We note that the ratio Rfs = ( 1 ^T 2 — l^-Pi)/(l^Pi — I^Pq) = 0.56(19) from our 
improved calculation is in better agreement with the experimental value (0.66). For the less 
accurate action we find 1.11(26) for this ratio. 

There will be additional uncertainties due to radiative corrections to the coefficients Cj 
in CNRQCD- Discretisation errors in the gluonic action will also effect our results and of 
course there will be corrections due to the introduction of sea quarks. Future studies into the 
accuracy of NRQGD will have to address these systematic errors. 
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Table 1: transformation behaviour of 66 -states in their non-relativistic notation 


State 


dA^mv^ 



0.4391(3) 

0.4416(3) 

0.00265(5) 

2^50 

0.688(23) 

0.679(22) 

- 

US! 

0.4497(4) 

0.4539(3) 

0.00439(6) 

DH 

0.687(23) 

0.681(20) 

- 

15 

0.4470(3) 

0.4508(3) 

0.00396(9) 

25 

0.69(3) 


- 


0.632(7) 

0.635(7) 

0.0033(1) 


0.619(7) 

0.621(7) 

0 . 0021 ( 2 ) 

3Pi 

0.628(7) 

0.627(8) 

0.0029(1) 

3pT2 

0.635(5) 

0.642(6) 

0.0053(1) 

^PE2 

0.632(7) 

0.642(6) 

0.0055(1) 

p 

0.630(3) 

0.634(4) 

0.0042(1) 


Table 2; Here we compare the absolute masses from different evolution equations and calculate 
their difference from a ratio fit with bootstrap errors. 
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Splitting 

0(mv'^) 

^^9 

0{mv‘^) 


Exp. value 

P - l^Si 

0.180(3) 

0.4398 

0.180(4) 

0.4398 

0.4398 GeV 


0.01073(5) 

0.0262(6) 

0.01227(11) 

0.0300(7) 

- 

23^1 - 135i 

0.24(2) 

0.59(5) 

0.23(2) 

0.56(5) 

0.5629 GeV 

liPi - 135i 

0.180(6) 

0.440(18) 

0.181(7) 

0.442(20) 

- 

13P2 - l^Po 

0.011(2) 

0.027(5) 

0.0143(34) 

0.035(8) 

0.0534 GeV 

13P2 - 

0.0045(14) 

0.011(3) 

0.0078(11) 

0.019(3) 

0.0213 GeV 

l3Pi - l3Po 

0.008(1) 

0.019(2) 

0.0070(13) 

0.017(3) 

0.0321 GeV 

13P2 - P 

0.00272(55) 

0.0066(14) 

0.0045(5) 

0.011(1) 

0.0131 GeV 

P - l^Pi 

0.0016(8) 

0.004(2) 

0.0036(6) 

0.0088(15) 

0.0082 GeV 

P - l^Po 

0.0078(11) 

0.019(3) 

0.0120(17) 

0.029(4) 

0.0403 GeV 

^kin 

3.84(8) 

9.38(25) 

3.97(8) 

9.70(29) 

9.4604 GeV 

Rsp 

1.33(11) 


1.28(11) 


1.2802 

Rfs 

0.56(19) 


1.11(26) 


0.6636 


Table 3; Summary of splittings. We show our results for both accuracies and convert them 
into physical units. The lattice spacing is determined from the IP — 1^5 splitting in both 
cases. RsP = {2^Si - l^Si)/{lP - Rfs = {^P 2 - ^Pi)l{^Pi - ^Po)- 


State 




E5o(p = (1,0,0)) 

0.458(1) 

0.462(1) 

0.00236(7) 

li5o(p = (1,1,0)) 

0.478(2) 

0.481(2) 

0.00212(8) 

l^*S'o(p = (1,1,1)) 

0.496(4) 

0.499(3) 

0.00191(10) 

li5o(p = (2,0,0)) 

0.522(4) 

0.524(5) 

0.0018(3) 

G>^kin 

3.77(7) 

3.93(7) 


T^5i(p = (1,0,0)) 

0.469(1) 

0.474(2) 

0.00427(5) 

i35i(p = (1 ,1,0)) 

0.489(2) 

0.493(2) 

0.00416(6) 

1^5'i(p = (1,1,1)) 

0.508(2) 

0.511(3) 

0.00406(7) 

i35i(p = (2 ,0,0)) 

0.531(5) 

0.535(5) 

0.00417(8) 

OiPIkin 

3.84(8) 

3.97(8) 



Table 4: The dispersion relation for ^Sq and ^Si. The momentum is given in lattice units of 
p = where L = 16. The last column is the difference for both evolution equations 

as obtained from a ratio fit. 
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Figure 1: Two examples of correlated fits for the ^5o-state, calculated with the improved 
action. The different symbols denote different choices for Ufu and the tmin dependence 
is shown. This is to demonstrate the consistency of our fit results as we change the fit 
prescription. We chose tmax = 48 for all fits. 
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Figure 2: Dispersion relation for and two different accuracies. The non-relativistic energy 
in lattice units is plotted vs. (||■)^|np where n = (l,0,0),n = (l,l,0),n = (1,1,1). From 
this we determine its kinetic mass to be 9.7(3) GeV (0(mu^)) and 9.4(3) GeV {0{mv^)) for 
arub = 1.71. 
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Upsilon Spectrum from NRQCD 
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Figure 3; Upsilon spectrum. We compare our results for two actions with different accura¬ 
cies. Crosses are the results from an action which is accurate to Circles denote an 

improved calculation where spin-corrections, C 4 , C 5 and cq, have been taken into account up to 
order 0{mv^). Bursts denote the experimental value where available. The splittings, relative 
to the are shown. The scale has been calculated from the IP — splitting. 



Figure 4: Upsilon fine structure. The splittings, relative to the P, are shown. Experimental 
values are shown as dashed lines. Crosses denote accuracy up to 0{mv^), circles up to 0{mv^). 
The scale has been calculated from the IP — splitting. 
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